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ABSTRACT 

Resistive relativist ic magnetohydrodynamic (RRMHD) simulations are applied to investigate the 
system evolution of relativistic magnetic reconnection. A time-split Harten-Lan-van Leer method is 
employed. Under a localized resistivity, the system exhibits a fast reconnection jet with an Alfvenic 
Lorentz factor inside a narrow Petschek-type exhaust. Various shock structures are resolved in and 
around the plasmoid such as the post-plasmoid vertical shocks and the "diamond-chain" structure 
due to multiple shock reflections. Under a uniform resistivity, Sweet-Parker- type reconnection slowly 
evolves. Under a current-dependent resistivity, plasmoids are repeatedly formed in an elongated 
current sheet. It is concluded that the resistivity model is of critical importance for RRMHD modeling 
of relativistic magnetic reconnection. 
Subject headings: magnetic reconnection — magnetohydrodynamics (MHD) — relativistic processes 
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Magnetic re connection (Sweet 1958 Parker 1957 
Pet schek|| 1964 ) is the driver ol explosive events m space, 
astrophysical, and laboratory plasmas. The reconnec- 
tion process attracts growing attentions to explain flar- 
ing events (Lyutikov 2006; Giannios et al. 2009) and the 



magn etic annihilation (Coroniti 1990; Lyubarsky & Kirk 
2001) in relativistic plasma environments. Basic prop- 
erties of relativistic magnetic reconnection have been 
discussed by relati vistic ma gnetohydrodynamic (R MHD) 
theories ( Blackman fc Fiel d 1994; Lyubarsky 2005 ; [Ten- 
Bar geet al ||2Ul(j| . liT^articular, |Lyubarsky| ( |20O51 rs 
careful work brought significant insights. 

In the last decade, modern simulation works have re- 
vealed many features of relativistic magnetic reconnec- 
tion. From the viewpoint of kinetic physics, it is widely 
recognized that reconnectio n is an efficient particle accel- 
erato r (Zenitani fc Hoshin^|2001 2007 Jaroschek et al. 
2004). It was further found that the system is cr ucially 



intluenced by the guide- field magnetic topology ( Zeni 
tani fc Ho shino 2005) and the radiative cooling ettects 



( [Jaroschek fc Hoshino]|2009 ). In a fluid scale, the recon 
nection syst em has been expl ored by a resistive RMHD 
(RRMHD) (Watanabe fc Yokoyama 2006) and relativis- 
tic tw o-fluid (Zenitani et al.„ 2009a b) models. ,ZenitanI| 
et al.| ([2009a) demonstrated and extensively analyzed a 
quasi-steady Petschek-type reconnection. 

There is a strong demand for further development of 
RRMHD reconnection work. Unlike the kinetic and two- 
fluid models, the RRMHD model is free from kinetic 
scales such as the skin depth and the gyro radius, and 
so it is highly desirable to study stellar-scale problems. 
In addition, plasmas are considered to be collision al in 



optically thick, radiati on-dominated environments (Uz- 
densky fc MacFadyen| |2006 ) . Since such plasmas can 
be approximated bva single RMHD fluid, the RRMHD 
code with non-FLEJ^ radiative transfer will be necessary 
for future modeling of radiative relativistic reconnection. 
However, no other RRMHD w ork has come out since 
Watanabe fc Yokoyamal (|2006|), because the RRMHD 



equations turned out to be numerically unstable. At- 
tempts to be tter solve th e RRMHD system have been 
undertaken. | Komissarov| (12007) pointed out that the 
non-ideal electric field behaves as stiff relaxation to the 
ideal RMHD condition. He split the equations into two 
parts: the stiff terms were analytically solved as an expo- 
nential decay, and the rest part was solved b y a two-step 
Harten-Lan-van Leer (HLL) method. Palenzuela et al. 
(2009) employed a hybrid scheme of an implicit Kunge- 
Kutta method for the e lectric field and an explicit one 



for the other variables. Dumbser fc Zanotti (2009) de- 



veloped a higher-order scheme in unstructured volumes. 
In this Letter, we present RRMHD simulations of rela- 
tivistic magnetic reconnection, based on recent advances 
in numerical schemes ( Komissarov^^2007} . 



2. NUMERICAL SETUP 

We employ the following RRMHD 
Lorentz-Heaviside notations with c = 1. 



equations in 



dtim^E xB)^V '{{p 



9^(7P) + V-(pn)=0 (1) 



-)I 



dt{£ 



B^^E^ 



^wuu - BB - Ee\ = (2) 

) + V-(m + £: X B) = (3) 



:0 



dtB^V xE-- 
dtE-V xB = -j 

j(e -^v X B — {E ■ v)v] =r]{j 



PcV) 



(4) 
(5) 
(6) 

(7) 



FLD is flux-limited diffusion 



In the above equations, p is the proper mass density, 
u = jv is the spatial part of the 4-velocity, m = jwu 
is the momentum, w is the enthalpy, £ = j'^w — p is 
the energy density, p is the proper pressure, pc is the 
charge density, and r] is the resistivity. We employ a F- 
law equation of state with an index of F = 4/3 and so 
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Table 1 

List of Simulation Runs 



Run Domain Grids 



At 
Ace 



CA, 























1 


120 X 120 


36002 


0.15 


0.1 


4.0 


0.894 


Eq. 


9 




IL 


180 X 180 


5400^ 


0.15 


0.1 


4.0 


0.894 


Eq. 


S 




lb 


120 X 120 


36002 


0.15 


0.1 


4.0 


0.894 


Eq. 


S 




Ic 


120 X 120 


36002 


0.15 


0.1 


4.0 


0.894 


Eq. 


TO 




2 


120 X 120 


24002 


0.2 


0.3 


1.33 


0.813 


Eq. 







3 


120 X 120 


24002 


0.2 


1.0 


0.4 


0.535 


Eq. 


§ 





Note. — Columns: (2) the domain size; (3) grid cells; (4) the 
timestep At/ Ax; (5) the background density pin] (6) the relevant 
magnetization ae^in] (7) Alfven speed CA,in', and (8) the resistivity 
model. 

the enthalpy w is given hy w = p -\- 4:p. In addition, we 
use a frame-independent magnetization parameter (j^ = 
b'^ /w = {B'^ — E'^)/w and the conductivity S = r]~^ in 
our discussion. S also stands for the magnetic Reynolds 
number based on the current sheet thickness L = 1 and 
the typical speed of c = 1, i.e., S = Lc/r]. 

We study two-dimensional system evolutions in the x- 
z plane. We employ a Harris-like model as an initial 
configuration: B = 5otanh(z) x, j = Bq cosh~'^ (z) y, 
p = pocosh~^{z) + Pin, P = P, Po = Bl/2 = 1, n = 0, 
and E = r]j. The background density (pin), the result- 
ing magnetization {(Js,in)^ and the Alfven speed (cA^in) 
are shown in Table IT] In particular, cr^^in is a measure 
of the system relativity, because the reconnection out- 
fiow is expected to be ~ CA,in = [crs,in/{^s,in + 1)]"^^^- 
The reference parame ters (run 1) are set similar to those 
in the earlier works (jWatanabe fc Yokoyama||2QQ6 



itani 
tEe 



et al. 



:s (iWa 

I. We 



2009a|). ^We set 

ter= (0,0). 



Zen- 



the reconnection pomt at 
A small vector potential 



origm 
SAy = 0.06^0 exp[—(x^ + ^^)/4] is imposed to trigger 
reconnection. Because of the symmetry of the recon- 
nection system, we consider point-symmetric conditions 
at X = 0: a property / satisfies /(O, z) = /(O, —z) or 
— /(O, —z). Neumann conditions of zero normal deriva- 
tives are employed at infiow {z = ±60 [±90 in run IL]) 
and outfiow {x = 120 [180]) boundaries. The normal 
magnetic fields {Bz) are adjusted to ensure V • S = at 

the infiow boundaries. 

We solv e the RRMHD equations by a variant of Komis- 



(2007)'s time-split HLL method. We solve a subset 



of the lull equations, because our specific configuration 
ensures By, E^, E^, Uy, pc, jx, 3z, and [E - v') = 0. 



Primitive va riables are analytically recovered ( [Zenitani 
et al. 2009a) and then interpolated by a second order 
monotonized central (MC) limiter. We use the hyper- 
bolic d ivergence cleaning m ethod for the solenoidal con- 
dition ( jDedner et aLl|2002D . 

We study three ditierent resistivity models — uniform, 
spatially localized, and current-dependent ones. 

7?(x,2;)=7?o±(^i 



vif^p)-- 



(8) 

7?0)cOsh-'(V^2+z2), (9) 

if < P'i'c) 






(10) 



Here, j^ 



J^J/i 



p^ is a frame-independent vari- 



able and Ic is a threshold. The baseline value is ^o = 
tJq^ = 80 in all three models. We mainly explore the 



second model, whose localized value is Si = rji = 10. 
The third model depends on the space-like current inten- 
sity, and the anomalous resistivity starts to work when 
the current is 5 times stronger than the Harris current, 
Ic = 5Bo/po. These resistivity models are handled by the 
analytic part of the numerical scheme. In the current- 
dependent case, t he partial Arnpere's law dfEy = —jy in 
the split part (see Komissarov ( [2007 ), section 3) has the 
following solution: 

(jyoexp{-So-ft) {\jy\<plc) .^^. 

1^2/0- ^(sgnj^o)p/c5'o7^ {\jy\>pQ^ ^ ^ 

where jyo is the initial state of j^y. 

3. RESULTS 

We study run 1 with the localized resistivity (Equation 
[9|. After an initial adjustment stage of t < 10, the re- 
connection process starts around the origin, and then it 
transfers the upstream magnetic fields faster and faster. 
The electric field Ey at the reconnection point (hereafter 
Ey), which stands for the fiux transfer speed of recon- 
nection, exceeds E* ^ O.IBq at t '^ 50. The normal- 
ized reconnection rate 7^ = E*/{cA,in'Bin') ^ Vin'/vout 
(using Vout ^ CA^in' and Ey^in' - E*) reaches its typ- 
ical value of 7^ '^ 0.145 at t ~ 100. Here, the up- 
stream properties with the subscript in' are evaluated at 
(x,z) = (0,20). Around this time, we recognize key sig- 
natures that are discussed later in this section, and main 
structures s tart to move at the nearly constant speeds 
( Ugai| 1995|). Therefore, we think that the reconnection 
is well developed after t > 100. 

Figure [l] shows the snapshots at t = 195. Panels 
(a),(b), and (c) show the out-of-plane current jy, the 
4- velocity Ux, and the magnetization a^, respectively. A 
fast reconnection jet travels outwar d inside a narro w ex- 
haust between a pair of slow shocks ( Petschek 1964 ) . The 
jet collides with a big magnetic island (the so-called plas- 
moid) in front of the current sheet further downstream. 
The typical outflow 4- velocity is ix^^ ~ 1.8 (7 ~ 2), and 
its maximum value is i^^^ = 2.75 (7 = 2.93) near the neck 
of the plasmoid {x ~ 80). The proper density inside the 
outflow exhaust is 3 times bigger than the that in up- 
stream region. The enthalpy flux {A^ypUx) carries 90% of 
outgoing energy flow at x = 60 due to the high pressure. 

Since the plasmoid suddenly compresses the surround- 
ing plasmas, the increased pressure pushes plasmas in 
the fleld-aligned directions. In the post-plasmoid regions 
(60 < X < 80), the plasma inflow is pushed in the —x- 
direction {ux < 0; blue regions in Figure [ijb)). In Fig- 
ures [ijb) and (c), a pair of vertical shocks can be seen 
at X ~ 60 between the left regions and the post-plasmoid 
regions. They are slow shocks and they move to the 
right at a constant speed of ~0.41c in the well-developed 
stage. Note that the left side is the shock downstream 
and the right is the upstream. In Figure [ijc), cr^ is en- 
hanced around the post-plasmoid regions. This is due 
to the lower pressure, caused by a fleld-aligned expan- 
sion of plasmas. The maximum speed in the outflow ex- 
haust is comparable with Alfven speed of the enhanced 
<j£ (~ 7.7) around there. Another pair of small shocks 
are found outside the plasmoid at x ~ 105, as indicated 
by the arrow in Figure flic) . 
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(a) Out-of-plane current (Jy) (b) Outflow 4-velocity (Ux=yVx) 
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Figure 1. Snapshots of run 1 at t = 195 in the x-z plane. Contour hues show the magnetic fields, (a) The out-of-plane 
(b) The 4-velocity Ux and the velocity vectors v. When \v\ = c, the length of the arrow is 6. (c) The magnetization de = 
(d) The 4-velocity Ux in 95 < x < 120 and — 4 < z < 4 (the small box in panel (b)). 



current jy/jo- 
{B^ - E^)/w. 



The plasmoid is surrounded by strong positive currents 
as can be seen in Figure [ij a). They are slow shocks (Ugai 



1995), which are connected to the Petschek slow shocks 



Due to the plasma heating and the magnetic energy dis- 
sipation across those slow shocks, the magnetization is 
weak inside the plasmoid (Figure [ijc)). One exception 
is the blue arc region, where the reconnected magnetic 
fields {Bz) are piled up at x '^ 84-89. There is a tan- 
gential discontinuity on the right vicinity of the piled-up 
region at x ^ 90. It looks like the small bright 'C- 
shaped structure in Figure fll a). There are a pair of re- 
verse current structures around x ~ 90-110 (blue regions 
in Figure [ija)). The magnetic fields are so bent across 
these oblique structures that their tangential components 
change the polarity. They are intermediate shocks (lAbei 
^nioshino 2001). 

"Figure flld) shows the i^^^-profile near the plasmoid 
edge. The same domain is indicated by the yellow box in 
Figure fljb). Plasmas are stationary around the central 
region (z ~ 0) because the high-density Harris sheet plas- 
mas are confined. Meanwhile, the twin plasmoid fiows 
travel fast in the upper and lower regions {z ~ ±3). They 
are separated by the intermediate shocks from the cen- 
tral region. The plasmoid edges propagate at ~0.75c and 
the plasma speed is ~0.66c there. Interestingly, we find 
the chain of diamond- shaped features around x ~ 102- 
115, z ~ 0. We call them the "diamond-chain" struc- 
ture. Since the twin plasmoid edges move faster than 
the fast mode or the sound speed inside the central region 
(~0.53c), bow shocks propagate from there, and then the 
shock fronts are refiected by the intermediate shocks mul- 
tiple times. The diamond-chain starts to develop in ear- 
lier stages of simulation and Figure [ijd) shows that six or 



more refiections occurred. We confirmed that the struc- 
ture is virtually unaffected by the boundary effect by 
carrying out a larger run (run IL). In the right side, the 
current sheet looks 60-100% wider than the initial con- 
dition. This is due to the current sheet diffusion whose 
timescale is L'^Sq ^ Sq = 80. 

Next, we study the dependence on the upstream mag- 
netization. We compare runs 1-3 in their well-developed 
stages. Typical reconnection rates are 71 = 0.145,0.125, 
and 0.102, respectively — reconnection evolves faster in 
higher-(j£ cases. Shown in Figure |2Va) are the 4- velocities 
in the outfiow exhaust. Since u^ anomalously increases 
near the plasmoid neck, we evaluate the typical 4-velocity 
in the left side of the vertical shocks. The typical 4- 
velocities (white squares) are in better agreement with 
the Alfven value y^ (dashed line) than the maximum 
(black squares). We also find that the upstream magne- 
tization gradually changes during the system evolution. 
For example, a^ decreased from 4 to 3.5 around x ^ 60 
in Figure mc). Considering this, the agreement becomes 
even better. This relation bridges the well-known nonrel- 
ativistic result of Vout ^ c^,in and the relativistic s caling 



in the hi^h-a^Ar^, lini it, ^out ^ crs,in ( [Lyubarsky 



2005 



Zenitani et al.||2009a[ ). 

Figure^^Jpresents the opening angle Opk of Petschek 
curren t layers. For reference, results of previous two-fiuid 
works ("Zenitani et al.i2009a') are presented too. In both 
cases, the angles become narrower as a^ increases. This 
intuitively fits the traditional explanation of the Petschek 
outfiow. In the jet frame, the outfiow exhaust expands 
by a local Alfven speed, (x B'^ ^ Bz/^out- In the ob- 
server's frame, such transverse motion in the relativis- 
tic jet looks slower by another factor of 7out- As a re- 
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suit, given that all the other conditions are the same, 
the Petschek angle is rescaled by ^~^^ ^ (1 + a^^^nj^'^- 
Assuniing that the slope angle of the upstream mag- 
netic field is Om ^ '7^, the Petschek angle would be 
a sizable fraction of 7^/(1 + cTe^in)- For an order esti- 
mate, the shaded region in Figure Mb) shows the range 
of 0.05 < {1 -\- a^^in) < 0.1. The angle analysis was 
further extended to ds^in = 9 and com pared with the 
high-a^^in theory in the two-fluid work (Zenitani et al. 
2009a|). To comprehensively discuss the outflow speed 



and the Petschek angle, a generalized theory which covers 
the full range from the nonrelativistic low-a^^in regime to 
the relativistic high-cr^^in limit needs to be developed. 



(a) Outflow 4-velocity 



(b) Opening angle 
0.1 




Sigma (a) 



Sigma (a) 



Figure 2. (a) Maximum and typical 4- velocities Ux as a function 
of the initial (Je,in- The dashed line indicates an Alfvenic value, 
^/o^. (b) Petschek openi ng angles (black line) and the relevant 
results in two-fluid work ( jZenitani et al.||2009a| (gray line). The 
shadow shows 0.05 < (1 -[- (Je^in) ^ < O.i. 



We further study the effect of the resistivity model. 
Run lb employs the uniform resistivity (Equation|8| with 
the same initial conditions as run 1. In this case, the 
system evolves very slowly. The plasmoid structure is 
visible after t ~ 400 and the outflow becomes relativistic 
Ux ^ I 8it t = 475-500. FigurepFa) shows the snapshot at 
t = 500. The reconnected current sheet remains thin and 
it is simply elongated in the x-direction: the system ex- 
hibits a Sweet-Parker-type reconnection. The reconnec- 
tion rate is slower by an order of magnitude, IZ ~ 0.01- 
0.02. Obviously the localized resistivity changed the sys- 
tem evolution. Run Ic employs the current-dependent 
resistivity (Equation 10). By this model, we intend to 
limit the current in low-dense regions, where the system 
does not have a lot of current carriers. The run starts 
from the intermediate data of run 1 at t = 25 to speed 
up the onset. In this case, a Sweet-Parker- type recon- 
nection similarly grows, but plasmoids repeatedly grow 
in a single elongated current sheet. Figure [sj^b) shows a 
well-developed stage at t = 215. One plasmoid stays near 
the origin because of the symmetric boundary condition. 
A new plasmoid starts to grow at x ~ 29 (indicated by 
the white arrow in Figure [sFb)), and the older one is 
ejected to the downstream, x ~ 55. Such a repeated for- 
mation of plasmoids is one of the ubiquitous features of 
reconnection. Their distances are comparable with the 
typical wavelength of the tearing mode ~ 20-30. The 
diamond-chain develops around the edge of the biggest 

t'asmoid, although it is difficult to distinguish in Figure 
b). The anomalous resistivity plays a role both in the 
Dugated current sheet and at the slow shocks surround- 
ing the post-plasmoid regions. The conductivity changes 



S ^ SO ^ 15-30 due to the intense currents there. Over- 
all evolution is slower than that of run 1, but faster than 
that of run lb. Using the electric field at the reconnec- 
tion point around x ^ 8, an equivalent reconnection rate 
is 7^ - 0.05-0.06. 

(a) Run lb (t=500) : Outflow 4-velocity (Ux=yVx) 

20 1 




20 

(b) Run Ic (t=215) 




Figure 3. Snapshots for (a) run lb at t = 500 and for (b) run 
Ic at t = 215: the 4-velocity Ux (color), the velocity v (arrows; the 
length is 6 when \v\ = c), and the magnetic fields (contour lines). 



4. DISCUSSION AND SUMMARY 
Our results corroborate the first RRMHD work by 



Watanabe & Yokoyama (2006) on many aspects such 
as Alfvenic outflow, the compression ratio, and the nar- 
rower slow shock structure. Judging from their inflow 
speed {vin/cA,in ^ ^)^ they observed faster reconnec- 
tion probably due to the stronger localized resistivity of 
Si ~ 0(1). Our resu lts also agree with the two-fluid 
work by Zenitani et al.| (2009a) such as the typical outflow 
parameters of the reference runs and the Petschek angle 
relations. Some visible differences are attributed to the 
two- fluid effects. The current layers look much sharper 
in this work, because the fluid inertial effects tend to 
smooth structures in the two-fluid model. The plasmoid 
looks longer in the x-direction. The outflow jet imme- 
diately hits the downstream plasmas in the x-direction 
in the RRMHD model, while plasma out-of-plane mo- 
tion softens the reflection in the two-fluid model. Con- 
cluding this and those works, (1) Alfvenic outflow, (2) 
increasing reconnection rate, and (3) narrower opening 
angle, can be regarded as common features of relativistic 
Petschek reconnection. Since a narrower outflow exhaust 
suppresses the energy throughput, (3) does not immedi- 
ately fit (2). However, Zenitani et al. (2009b) pointed out 
that the relativistic enthalpy tlux transports huge energy 
per a cross section even in a narrower exhaust. 
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Thanks to the stable RRMHD code, we successfuhy 
resolved shock structures in and around the plasmoid. 
To our knowledge, the post-plasmoid vertical shocks and 
the diamond-chain structure are new discoveries. The 
diamond chain develops when the plasmoid edge speed 
(a sizable fraction of CA^in) exceeds the sound speed in 
the neutral Harris sheet. We predict that it also devel- 
ops in nonrelativistic MHD simulations once the above 
condition is met. Since multiple shocks are confined in a 
narrow region, the plasmoid edge can be a potential site 
of shock acceleration of particles. 

From the numerical viewpoint, relativistic reconnec- 
tion in the high-a^ regimes is challenging for an explicit- 
type RRMHD code. Since the analytic part of the 
scheme contains the time constant of —Sj (e.g.. Equa- 
tion [TT^), employing Alfvenic outfiow speed, a code re- 
quires At < 'S'o""'^(l + crs^in)~^^'^' As shown in Section 
3, the outfiow is faster near the plasmoid neck, and 
so the restriction becomes even severer. Thus, high 
time (and appropriate spatial) resolutions are necessary 
in high-(j£ and moderate-^ regimes of our interest. In 
addition, a plasma solution is not always stable even 
when the total momentum and energy are well handled. 
Let us extract the plasma momentum from Equation [2J 
dtm-\-\/ ■ {pi + wuu) = j X B. Focusing on the source 
term which is usually problematic, we see that the plasma 
stability condition, \Am\/w -^ \j x B\At/w < (9(0.1), 
will be critical in higher-a^ regimes, unless Ohm's law 
strongly limits the current. The se issues need to be ini - 
proved by an implicit approach ( Palenzuela et al.||2QQ9 ). 

From the physics viewpoint, those restrictions obvi- 
ously come from an immature form of Ohm's law. Note 
that the RRMHD model contains the displacement cur- 
rent dfE in Ampere's law (Equation Isl). Since it oper- 
ates in a short timescale of the plasma frequency, the 
nonrelativistic MHD code usually drops it. However, the 
RRMHD system needs it for the energy and momentum 
conservation. On the other hand, the RRMHD uses a 
simple Ohm's law with the scalar resistivity (Equation[7l. 
This is a time-stationary form (Blackman fc Fie ld 1953) 
and there is no general consensus on extended forms. 
Combining the two equations with different timescales, 
the RRMHD system inevitably becomes problematic, es- 
pecially in its impulsive phase. Thus , a f urther devel- 
opment of Ohm's law (e.g., Gedahnl ( |1996 )) or the re- 
construction of the entire RRMHD system (e.g., iKoidel 
( [2009) ) is highly desirable. 

In order to study large-scale problems which contain 
reconnect ion in arbitrary locations, we can no longer use 
the localized resistivity at a fixed location. We do need 



good physical or phenomenological anomalous resistivity 
models. Our attempt is a first step toward this direction, 
and the current-dependent resistivity model did exhibit 
different system evolution. This clearly tell us that the 
modeling of the effective resistivity, a long-standing prob- 
lem of the entire reconnection research, is very important 
for RRMHD modeling of relativistic magnetic reconnec- 
tion as well. 

In summary, we explored an RRMHD model of rel- 
ativistic magnetic reconnection. We corroborated ear- 
lier works and further studied fine plasmoid structures. 
Three different resistivity models are examined and they 
showed different system evolutions. It is crucially impor- 
tant to model an effective resistivity for RRMHD mod- 
eling of magnetic reconnection. 
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